## Loading packages ###
packages <- c("haven", "multiwayvcov", "lmtest", 
              "MatchIt", "dplyr", "cobalt", "ggpubr",  "separationplot", "interplot",
              "texreg", "MASS", "konfound", "stargazer")

sapply(packages, function(x){library(x, character.only = TRUE)})

## Loading data ##
df <- read_dta("Data/200731_data file participation compliance CLRTAPwithIGOcount.dta")


## Defining marginal effect plot function. Courtesy of Miles D. Williams: https://rpubs.com/milesdwilliams15/381372
meplot <- function(model,var1,var2,int,vcov,ci=.95,
                   xlab=var2,ylab=paste("Marginal Effect of",var1),
                   main="Marginal Effect Plot",
                   me_lty=1,me_lwd=1,me_col="black",
                   ci_lty=1,ci_lwd=.5,ci_col="black",
                   yint_lty=2,yint_lwd=1,yint_col="black"){
  alpha <- 1-ci
  z <- qnorm(1-alpha/2)
  beta.hat <- coef(model)
  cov <- vcov
  z0 <- seq(min(model.frame(model)[,var2],na.rm=T),max(model.frame(model)[,var2],na.rm=T),length.out=1000)
  dy.dx <- beta.hat[var1] + beta.hat[int]*z0
  se.dy.dx <- sqrt(cov[var1,var1] + z0^2*cov[nrow(cov),ncol(cov)] + 2*z0*cov[var1,ncol(cov)])
  upr <- dy.dx + z*se.dy.dx
  lwr <- dy.dx - z*se.dy.dx
  ggplot(data=NULL,aes(x=z0, y=dy.dx)) +
    labs(x=xlab,y=ylab,title=main) +
    geom_line(aes(z0, dy.dx),size = me_lwd, 
              linetype = me_lty, 
              color = me_col) +
    geom_line(aes(z0, lwr), size = ci_lwd, 
              linetype = ci_lty, 
              color = ci_col) +
    geom_line(aes(z0, upr), size = ci_lwd, 
              linetype = ci_lty, 
              color = ci_col) +
    geom_hline(yintercept=0,linetype=yint_lty,
               size=yint_lwd,
               color=yint_col)
}

#### Logit models ####
logit1 <- glm(Compl_dich ~ Participant, family = binomial(link = "logit"), data = df)
logit1$robust <- coeftest(logit1, vcov. =  cluster.vcov(logit1, df$Country))


logit2 <- glm(Compl_dich ~ Participant + 
                Ambition_level, family = binomial(link = "logit"), data = df)
logit2$robust <- coeftest(logit2, vcov. = cluster.vcov(logit2, df$Country))

logit3 <- glm(Compl_dich ~ 
                Participant +
                Ambition_level +
                Govt_Eff_WGI, family = binomial(link = "logit"), data = df)
logit3$robust <- coeftest(logit3, vcov. = cluster.vcov(logit3, df$Country))

logit4 <- glm(Compl_dich ~ Participant + 
                Ambition_level +
                Govt_Eff_WGI +
                East_Europe
                , family = binomial(link = "logit"), data = df)
logit4$robust <- coeftest(logit4, vcov. = cluster.vcov(logit4, df$Country))


logit5 <- glm(Compl_dich ~ Participant + 
                Ambition_level +
                Govt_Eff_WGI +
                East_Europe +
                EU_member
              , family = binomial(link = "logit"), data = df)
logit5$robust <- coeftest(logit5, vcov. = cluster.vcov(logit5, df$Country))

logit6 <- glm(Compl_dich ~ Participant + 
                Ambition_level +
                Govt_Eff_WGI +
                East_Europe +
                EU_member +
                GDP_growth
              , family = binomial(link = "logit"), data = df)
logit6$robust <- coeftest(logit6, vcov. = cluster.vcov(logit6, df$Country))

logit7 <- glm(Compl_dich ~ Participant + 
                Ambition_level +
                Govt_Eff_WGI +
                East_Europe +
                EU_member +
                GDP_growth+
                govt_autonomy
              , family = binomial(link = "logit"), data = df)
logit7$robust <- coeftest(logit7, vcov. = cluster.vcov(logit7, df$Country))


logit8 <- glm(Compl_dich ~ Participant *
                Ambition_level + 
                Govt_Eff_WGI + 
                East_Europe + 
                EU_member +
                GDP_growth+
                govt_autonomy,
              family = binomial(link = "logit"), data = df)
logit8$robust <- coeftest(logit8, vcov. = cluster.vcov(logit8, df$Country))




stargazer(logit1, logit2, logit3, logit4, logit5, logit6, logit7, logit8, 
          se = list(logit1$robust[,2], 
                    logit2$robust[,2], 
                    logit3$robust[,2], 
                    logit4$robust[,2], 
                    logit5$robust[,2], 
                    logit6$robust[,2], 
                    logit7$robust[,2], 
                    logit8$robust[,2] ), 
          p = list(logit1$robust[,4], 
                   logit2$robust[,4], 
                   logit3$robust[,4], 
                   logit4$robust[,4], 
                   logit5$robust[,4], 
                   logit6$robust[,4], 
                   logit7$robust[,4], 
                   logit8$robust[,4] ), 
          type = "html", 
          title = "Logit regression coefficients (standard errors). Dichotomous DV", 
          covariate.labels = c("Participation",
                               "Ambition",
                               "Government effectiveness",
                               "Eastern Europe",
                               "EU member",
                               "GDP growth",
                               "Veto players", 
                               "Participation*Ambition"),
          dep.var.caption = "", 
          dep.var.labels.include = FALSE, 
          notes = "Standard errors clustered on countries.", 
          notes.align = "l", 
          notes.append = TRUE, 
          out = "logitTable.rtf", 
          font.size = "scriptsize")


meplot(logit8, 
       var1 = "Participant", 
       var2 = "Ambition_level", 
       int="Participant:Ambition_level", 
       vcov = cluster.vcov(logit8, df$Country), 
       main = "")+
  xlab("Ambition")+
  ylab("Marginal effect of Participation")+
  theme_classic()
ggsave("InteractionLogit8.png",  height = 5, width = 5)







set.seed(4556)
simb <- mvrnorm(n = 1000, 
                     mu = coefficients(logit7),
                     Sigma = cluster.vcov(logit7, df$Country))

setX <- cbind(1,
              c(0,1),
              mean(df$Ambition_level, na.rm = TRUE), 
              mean(df$Govt_Eff_WGI, na.rm = TRUE), 
              0, #East_Europe
              1, #EU_member
              mean(df$GDP_growth, na.rm = TRUE),
              mean(df$govt_autonomy, na.rm = TRUE))
xBeta <- setX %*% t(simb)
expY <- 1/(1+exp(-xBeta))
expY <- as.data.frame(t(expY))
colnames(expY) <- c("Non-Participant", "Participant")
expY$Difference <-  expY$`Participant` - expY$`Non-Participant`


simulated <- as.data.frame(cbind(c("Nonparticipant", "Participant", "Difference"), 
                   t(apply(X = expY, MARGIN = 2, FUN = quantile, probs = c(.025,.5,.975)))))
names(simulated) <- c("Participant", "lower", "point", "upper")
simulated <- mutate(simulated, across(.cols = c(lower,point,upper), .fns = as.numeric))

ggplot(dplyr::filter(simulated, Participant != "Difference"), 
       aes(x = Participant, y = point, ymin = lower, ymax = upper))+
  geom_errorbar(width = 0.1) +
  geom_point()+
  theme_classic()+
  ylim(0.3, 1)+
  xlab("")+
  ylab("Predicted probability of compliance")
ggsave("predictionsLogit7.png",  height = 5, width = 5)



ols1 <- lm(Compl_cont ~ Participant, data = df)
ols1$robust <- coeftest(ols1, vcov. = cluster.vcov(ols1, df$Country))

ols2 <- lm(Compl_cont ~ Participant + 
             Ambition_level, data = df)
ols2$robust <- coeftest(ols2, vcov. = cluster.vcov(ols2, df$Country))

ols3 <- lm(Compl_cont ~ Participant + 
             Ambition_level+
             Govt_Eff_WGI, data = df)
ols3$robust <- coeftest(ols3, vcov. = cluster.vcov(ols3, df$Country))

ols4 <- lm(Compl_cont ~ Participant + 
             Ambition_level+
             Govt_Eff_WGI+
             East_Europe, data = df)
ols4$robust <- coeftest(ols4, vcov. = cluster.vcov(ols4, df$Country))


ols5 <- lm(Compl_cont ~ Participant + 
             Ambition_level+
             Govt_Eff_WGI+
             East_Europe+
             EU_member, data = df)
ols5$robust <- coeftest(ols5, vcov. = cluster.vcov(ols5, df$Country))

ols6 <- lm(Compl_cont ~ Participant + 
             Ambition_level+
             Govt_Eff_WGI+
             East_Europe+
             EU_member+
             GDP_growth, data = df)
ols6$robust <- coeftest(ols6, vcov. = cluster.vcov(ols6, df$Country))

ols7 <- lm(Compl_cont ~ Participant + 
             Ambition_level+
             Govt_Eff_WGI+
             East_Europe+
             EU_member+
             GDP_growth+
             govt_autonomy, data = df)
ols7$robust <- coeftest(ols7, vcov. = cluster.vcov(ols7, df$Country))


ols8 <- lm(Compl_cont ~ Participant * 
             Ambition_level+
             Govt_Eff_WGI+
             East_Europe+
             EU_member+
             GDP_growth+
             govt_autonomy, data = df)
ols8$robust <- coeftest(ols8, vcov. = cluster.vcov(ols8, df$Country))


stargazer(ols1,ols2,  ols3, 
        ols4, ols5, ols6, 
         ols7,  ols8,
          se = list(ols1$robust[,2], 
                    ols2$robust[,2], 
                    ols3$robust[,2], 
                    ols4$robust[,2], 
                    ols5$robust[,2], 
                    ols6$robust[,2], 
                    ols7$robust[,2], 
                    ols8$robust[,2] ), 
          p = list(
                   ols1$robust[,4], 
                   ols2$robust[,4], 
                   ols3$robust[,4], 
                   ols4$robust[,4], 
                   ols5$robust[,4], 
                   ols6$robust[,4], 
                   ols7$robust[,4], 
                   ols8$robust[,4]), 
          type = "html", 
          title = "Table 4 OLS regression coefficients (standard errors). Continuous DV.", 
          covariate.labels = c("Participation", 
                               "Ambition", 
                               "Government effectiveness", 
                               "Eastern Europe", 
                               "EU member", 
                               "GDP growth", 
                               "Veto players", 
                               "Participation * Ambition"), 
          dep.var.caption = "", 
          dep.var.labels.include = FALSE, 
          notes = "Standard errors clustered on countries.", 
          notes.align = "l", 
          notes.append = TRUE, 
          keep.stat = c("n", "ll", "aic", "rsq"),
          out = "OlsTable.rtf")



meplot(ols8, 
       var1 = "Participant", 
       var2 = "Ambition_level", 
       int="Participant:Ambition_level", 
       vcov = cluster.vcov(ols8, df$Country), 
       main = "")+
  xlab("Ambition")+
  ylab("Marginal effect of Participation")+
  theme_classic()
ggsave("InteractionOls8.png",  height = 5, width = 5)








###### Fixed effects. 

logit9 <-  glm(Compl_dich ~ Participant +
                 Ambition_level + 
                 Govt_Eff_WGI + 
                 East_Europe + 
                 EU_member +
                 GDP_growth +
                 govt_autonomy+ 
                 Oslo_prot +
                 Oslo_prot_05target +
                 Gothenburg_prot,
               family = binomial(link = "logit"), 
               data = df)

logit9$robust <- coeftest(logit9, cluster.vcov(logit9, df$Country))

logit10 <-  glm(Compl_dich ~ Participant +
                 Ambition_level + 
                 Govt_Eff_WGI + 
                 East_Europe + 
                 EU_member +
                 GDP_growth +
                 govt_autonomy+ 
                 Oslo_prot +
                 Oslo_prot_05target +
                 Gothenburg_prot +
                  as.factor(Country),
               family = quasibinomial(link = "logit"), 
               data = df)

logit10$robust <- coeftest(logit10, cluster.vcov(logit10, df$Country))


ols9 <- lm(Compl_cont ~ Participant +
             Ambition_level + 
             Govt_Eff_WGI + 
             East_Europe + 
             EU_member +
             GDP_growth +
             govt_autonomy+
             Oslo_prot +
             Oslo_prot_05target +
             Gothenburg_prot,
           data = df)
ols9$robust <- coeftest(ols9, cluster.vcov(ols9, df$Country))
ols10 <- lm(Compl_cont ~ Participant +
             Ambition_level + 
             Govt_Eff_WGI + 
             East_Europe + 
             EU_member +
             GDP_growth +
              govt_autonomy+
             Oslo_prot +
             Oslo_prot_05target +
             Gothenburg_prot +
             as.factor(Country),
           data = df)

ols10$robust <- coeftest(ols10, cluster.vcov(ols10, df$Country))


stargazer(logit9,
          ols9, ols10, 
          out = "ModelsFixedEffects.rtf",
          se  = list(logit9$robust[,2],
                     ols9$robust[,2], 
                     ols10$robust[,2]), 
          p = list(logit9$robust[,4],
                   ols9$robust[,4], 
                   ols10$robust[,4]), 
          type = "html", 
          title = "Fixed-effects logistic regression and OLS coefficients (standard errors)", 
          covariate.labels = c("Participation", 
                               "Ambition", 
                               "Government effectiveness", 
                               "Eastern Europe", 
                               "EU member", 
                               "GDP growth", 
                               "Veto players"), 
          dep.var.caption = "", 
          dep.var.labels.include = FALSE, 
          omit = c("_prot", "factor"), 
          omit.labels = c("Protocol fixed effects", "Country fixed effects"),
          notes = "Standard errors clustered on countries.", 
          notes.align = "l", 
          notes.append = TRUE, 
          keep.stat = c("n", "ll", "aic", "rsq"))

#### Matching ###

matchingData <- df %>% 
  dplyr::select(Compl_cont, 
                Compl_dich, 
                Participant,
                Ambition_level,
                Govt_Eff_WGI,
                East_Europe, 
                EU_member,
                GDP_growth,
                govt_autonomy,
                Helsinki_prot,
                Oslo_prot,
                Oslo_prot_05target,
                Gothenburg_prot,
                Country) %>% 
  na.omit()
mahalanobisFull <- matchit(Participant ~
                                  Ambition_level + 
                                  Govt_Eff_WGI + 
                                  East_Europe + 
                                  EU_member +
                                  GDP_growth +
                                  govt_autonomy +
                                  # Helsinki_prot +
                                  Oslo_prot +
                                  Oslo_prot_05target +
                                  Gothenburg_prot,
                                method = "nearest", 
                                distance = "mahalanobis",
                                replace = TRUE, 
                                data = matchingDataNoDep)
mahalanobisBalance <- love.plot(mahalanobisFull, stats = "mean.diffs", abs = FALSE, 
                                     var.names = c(Ambition_level = "Ambition level", 
                                                   Govt_Eff_WGI = "Government effectiveness", 
                                                   East_Europe = "Eastern Europe", 
                                                   EU_member = "EU member", 
                                                   GDP_growth = "GDP growth", 
                                                   govt_autonomy = "Political constraints", 

                                                   Oslo_prot = "Oslo Protocol", 
                                                   Oslo_prot_05target = "Oslo Protocol 2005 target", 
                                                   Gothenburg_prot = "Gothenburg Protocol"), 
                                     stars = "std", shapes = c(19, 17))+
  labs(subtitle = "Excluding domestic dep.", title ="")+
  xlim(-2.1, 1.3)

ggsave("mahalanobisBalance.png",  height = 10, width = 10)

matchedDataMahalanobis <- match.data(mahalanobisFull)

logitFullmatchedmahalanobis <- glm(Compl_dich ~
                                          Participant +
                                          Ambition_level +
                                          Govt_Eff_WGI +
                                          East_Europe +
                                          EU_member +
                                          GDP_growth +
                                          govt_autonomy +
                                          Oslo_prot +
                                          Oslo_prot_05target +
                                          Gothenburg_prot,
                                        family = binomial(link = "logit"),
                                        data = matchedDataMahalanobis, 
                                        weights = matchedDataMahalanobis$weight)
logitFullmatchedmahalanobis$robust <- coeftest(logitFullmatchedmahalanobis, vcov. = cluster.vcov(logitFullmatchedmahalanobis, matchedDataMahalanobis$Country))



olsFullmatchedmahalanobis <- lm(Compl_cont ~
                                       Participant +
                                       Ambition_level +
                                       Govt_Eff_WGI +
                                       East_Europe +
                                       EU_member +
                                       GDP_growth +
                                       govt_autonomy +
                                       Oslo_prot +
                                       Oslo_prot_05target +
                                       Gothenburg_prot,
                                     data = matchedDataMahalanobis, weights = matchedDataMahalanobis$weight)
olsFullmatchedmahalanobis$robust <- coeftest(olsFullmatchedmahalanobis, vcov. = cluster.vcov(olsFullmatchedmahalanobis, matchedDataMahalanobis$Country))




stargazer(logitFullmatchedmahalanobis, 
          olsFullmatchedmahalanobis, 
          out = "matchedModels.rtf",
          se  = list(logitFullmatchedmahalanobis$robust[,2],
                     olsFullmatchedmahalanobis$robust[,2]), 
          p = list(logitFullmatchedmahalanobis$robust[,4], 
                   olsFullmatchedmahalanobis$robust[,4]),  
          type = "html", 
          title = "Models estimated after matching. Logistic regression and OLS coefficients (standard errors)", 
          covariate.labels = c("Participation", 
                               "Ambition", 
                               "Government effectiveness", 
                               "Eastern Europe", 
                               "EU member", 
                               "GDP growth", 
                               "Veto players"), 
          dep.var.caption = "", 
          dep.var.labels.include = FALSE, 
          omit = c("_prot"), 
          omit.labels = c("Protocol fixed effects"),
          notes = "Standard errors clustered on countries.", 
          notes.align = "l", 
          notes.append = TRUE, 
          keep.stat = c("n", "ll", "aic", "rsq"))




